Preparation

load packages

library(brms)

## Loading required package: Rcpp

## Loading 'brms' package (version 2.21.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').

## 
## Attaching package: 'brms'

## The following object is masked from 'package:stats':
## 
##     ar

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(sjPlot)
library(performance)
library(rstan)

## Loading required package: StanHeaders
## 
## rstan version 2.32.6 (Stan version 2.32.2)
## 
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)
## 
## 
## Attaching package: 'rstan'
## 
## The following object is masked from 'package:tidyr':
## 
##     extract

library(bayesplot)

## This is bayesplot version 1.11.1
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
##    * Does _not_ affect other ggplot2 plots
##    * See ?bayesplot_theme_set for details on theme setting
## 
## Attaching package: 'bayesplot'
## 
## The following object is masked from 'package:brms':
## 
##     rhat

library(loo)

## This is loo version 2.8.0
## - Online documentation and vignettes at mc-stan.org/loo
## - As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session. 
## 
## Attaching package: 'loo'
## 
## The following object is masked from 'package:rstan':
## 
##     loo

library(bayestestR)
library(psych)

## 
## Attaching package: 'psych'
## 
## The following object is masked from 'package:rstan':
## 
##     lookup
## 
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
## 
## The following object is masked from 'package:brms':
## 
##     cs

library(reshape2)

## 
## Attaching package: 'reshape2'
## 
## The following object is masked from 'package:tidyr':
## 
##     smiths

library(GGally)

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

library(ggridges)
library(DT)
library(knitr)
library(kableExtra)

## 
## Attaching package: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows

sessionInfo()

## R version 4.4.0 (2024-04-24)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sonoma 14.6.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Asia/Tokyo
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] kableExtra_1.4.0    knitr_1.48          DT_0.33            
##  [4] ggridges_0.5.6      GGally_2.2.1        reshape2_1.4.4     
##  [7] psych_2.4.6.26      bayestestR_0.14.0   loo_2.8.0          
## [10] bayesplot_1.11.1    rstan_2.32.6        StanHeaders_2.32.10
## [13] performance_0.12.2  sjPlot_2.8.16       lubridate_1.9.3    
## [16] forcats_1.0.0       stringr_1.5.1       dplyr_1.1.4        
## [19] purrr_1.0.2         readr_2.1.5         tidyr_1.3.1        
## [22] tibble_3.2.1        ggplot2_3.5.1       tidyverse_2.0.0    
## [25] brms_2.21.0         Rcpp_1.0.13        
## 
## loaded via a namespace (and not attached):
##  [1] mnormt_2.1.1         gridExtra_2.3        inline_0.3.19       
##  [4] sandwich_3.1-0       rlang_1.1.4          magrittr_2.0.3      
##  [7] multcomp_1.4-26      matrixStats_1.3.0    compiler_4.4.0      
## [10] systemfonts_1.1.0    vctrs_0.6.5          pkgconfig_2.0.3     
## [13] fastmap_1.2.0        backports_1.5.0      rmdformats_1.0.4    
## [16] utf8_1.2.4           rmarkdown_2.28       tzdb_0.4.0          
## [19] xfun_0.47            cachem_1.1.0         jsonlite_1.8.8      
## [22] sjmisc_2.8.10        ggeffects_1.7.0      parallel_4.4.0      
## [25] R6_2.5.1             bslib_0.8.0          stringi_1.8.4       
## [28] RColorBrewer_1.1-3   jquerylib_0.1.4      estimability_1.5.1  
## [31] bookdown_0.40        zoo_1.8-12           Matrix_1.7-0        
## [34] splines_4.4.0        timechange_0.3.0     tidyselect_1.2.1    
## [37] rstudioapi_0.16.0    abind_1.4-5          yaml_2.3.10         
## [40] codetools_0.2-20     sjlabelled_1.2.0     curl_5.2.1          
## [43] pkgbuild_1.4.4       lattice_0.22-6       plyr_1.8.9          
## [46] withr_3.0.1          bridgesampling_1.1-2 posterior_1.6.0     
## [49] coda_0.19-4.1        evaluate_0.24.0      survival_3.7-0      
## [52] ggstats_0.6.0        RcppParallel_5.1.9   xml2_1.3.6          
## [55] pillar_1.9.0         tensorA_0.36.2.1     checkmate_2.3.2     
## [58] stats4_4.4.0         insight_0.20.3       distributional_0.4.0
## [61] generics_0.1.3       hms_1.1.3            rstantools_2.4.0    
## [64] munsell_0.5.1        scales_1.3.0         xtable_1.8-4        
## [67] glue_1.7.0           emmeans_1.10.3       tools_4.4.0         
## [70] mvtnorm_1.2-6        grid_4.4.0           QuickJSR_1.3.1      
## [73] datawizard_0.12.2    colorspace_2.1-1     nlme_3.1-166        
## [76] cli_3.6.3            fansi_1.0.6          viridisLite_0.4.2   
## [79] svglite_2.1.3        Brobdingnag_1.2-9    sjstats_0.19.0      
## [82] V8_5.0.0             gtable_0.3.5         sass_0.4.9          
## [85] digest_0.6.37        TH.data_1.1-2        htmlwidgets_1.6.4   
## [88] htmltools_0.5.8.1    lifecycle_1.0.4      MASS_7.3-61

Set the theme

theme_set(theme_bw())
set.seed(123)
rstan_options(auto_write = T)
options(mc.cores=parallel::detectCores())

load data

dat <- read_csv("YesNo.csv")

## Rows: 99 Columns: 181
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (180): accumulation, accurate, adaptation, adjacent, allocation, alter, ...
## dbl   (1): Participant
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

convert to long format

data_long <- dat %>% 
  pivot_longer(
    cols = -Participant,
    names_to = "Item",
    values_to = "Response") %>%
  mutate(
    Participant = factor(Participant))

make correct answers dataset

item_list <- c("accumulation",  "accurate", "adaptation",   "adjacent", "allocation",   "alter",    "amendment",    "arbitrary",    "assembly", "assurance",    "aware",    "behalf",   "capable",  "clarity",  "coincide", "colleagues",   "comprise", "concurrent",   "consultation", "contradiction",    "conversely",   "crucial",  "denote",   "depression",   "deviation",    "discretion",   "discrimination",   "diversity",    "domain",   "duration", "edition",  "eliminate",    "enable",   "energy",   "enforcement",  "enormous", "equipment",    "equivalent",   "erosion",  "estate",   "ethical",  "exceed",   "explicit", "exploitation", "exposure", "external", "flexibility",  "forthcoming",  "foundation",   "gender",   "guidelines",   "hierarchical", "incentive",    "incidence",    "inclination",  "incompatible", "inherent", "inhibition",   "initiatives",  "innovation",   "insights", "inspection",   "instructions", "integrity",    "intelligence", "intensity",    "intervention", "intrinsic",    "invoke",   "likewise", "logic",    "manipulation", "mediation",    "minimal",  "mode", "nonetheless",  "nuclear",  "odd",  "ongoing",  "orientation",  "persistent",   "practitioners",    "protocol", "publication",  "pursue",   "quotation",    "reluctant",    "restore",  "scenario", "simulation",   "sphere",   "stability",    "straightforward",  "submitted",    "substitution", "supplementary",    "symbolic", "thereby",  "thesis",   "topic",    "transmission", "undergo",  "unique",   "vehicle",  "visible",  "vision",   "visual",   "welfare",  "absolvention", "ackery",   "almanical",    "amagran",  "amphlett", "annobile", "ashment",  "asslam",   "atribus",  "attard",   "bastionate",   "benevolate",   "berrow",   "cambule",  "captivise",    "carpin",   "causticate",   "charactal",    "coath",    "coppard",  "crucialate",   "cymballic",    "decorite", "defunctionary",    "deliction",    "descript", "doole",    "dring",    "dyment",   "ebullible",    "eluctant", "expostulant",  "galpin",   "hapgood",  "hawther",  "hegedoxy", "hermantic",    "hignall",  "hislop",   "hoult",    "interisation", "jemmett",  "joice",    "keable",   "kearle",   "loveridge",    "majury",   "mastaphitis",  "neutration",   "nichee",   "oestrogeny",   "paralogue",    "perceptacle",  "peritonic",    "prowt",    "rainish",  "remonic",  "samphirate",   "savery",   "savourite",    "stattock", "shuddery", "tamnifest",    "tearle",   "tindle",   "transcendiary",    "vardy",    "viggers",  "warman",   "waygood",  "whitrow",  "wintle")
correct_answer_list <- c("Yes", "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "Yes",  "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No",   "No")

correct_answers <- data.frame(
  Item = item_list,
  Correct_Answer = correct_answer_list)

combine data

combined_data <- merge(data_long, correct_answers, by = "Item")

Scoring Methods

• The Hit column is 1 when Response is “Yes” and Correct_Answer is “Yes”

• The Miss column is 1 when Response is “No” and Correct_Answer is “Yes”

• The False_Alarm column is 1 when Response is “Yes” and Correct_Answer is “No”

• The Correct_Rejection column is 1 when Response is “No” and Correct_Answer is “No”

• Scoring_1: 1 point if Hit is 1, 0 points for Correct_Rejection, Miss, and False_Alarm

• Scoring_2: 1 point if Hit is 1, 1 point if Correct_Rejection is 1, 0 points for Miss and False_Alarm

• Scoring_3: 2 points if Hit is 1, 1 point if Correct_Rejection is 1, 0 points for Miss and False_Alarm

• Scoring_4: 3 points if Hit is 1, 2 points if Correct_Rejection is 1, 1 point if Miss is 1, 0 points if False_Alarm is 1

scoring_data <- combined_data %>%
  mutate(
    Hit = ifelse(Response == "Yes" & Correct_Answer == "Yes", 1, 0),
    Miss = ifelse(Response == "No" & Correct_Answer == "Yes", 1, 0),
    False_Alarm = ifelse(Response == "Yes" & Correct_Answer == "No", 1, 0),
    Correct_Rejection = ifelse(Response == "No" & Correct_Answer == "No", 1, 0),
    Scoring_1 = ifelse(Hit == 1, 1, 0),
    Scoring_2 = ifelse(Hit == 1 | Correct_Rejection == 1, 1, 0),
    Scoring_3 = ifelse(Hit == 1, 2, ifelse(Correct_Rejection == 1, 1, 0)),
    Scoring_4 = ifelse(Hit == 1, 2, ifelse(Correct_Rejection == 1, 1, ifelse(False_Alarm == 1, 0, 1))),
    Scoring_5 = ifelse(Hit == 1, 3, ifelse(Correct_Rejection == 1, 2, ifelse(False_Alarm == 1, 0, 1))))

convert to long format

long_format_data <- scoring_data %>%
  pivot_longer(
    cols = c(Scoring_1, Scoring_2, Scoring_3, Scoring_4, Scoring_5),
    names_to = "Scoring_Method",
    values_to = "Score")

long_format_data <- mutate(
  long_format_data,
  Participant = factor(Participant),
  Scoring_Method = factor(Scoring_Method))

Descriptive Statistics

participant_scores <- long_format_data %>%
  group_by(Participant, Scoring_Method) %>%
  summarise(Total_Score = sum(Score, na.rm = T)) %>%
  pivot_wider(names_from = Scoring_Method, values_from = Total_Score)

## `summarise()` has grouped output by 'Participant'. You can override using the
## `.groups` argument.

participant_scores %>%
  dplyr::select(Scoring_1, Scoring_2, Scoring_3, Scoring_4, Scoring_5) %>% 
  describe() %>% 
  kable(format = "html") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = F)

## Adding missing grouping variables: `Participant`

long_data <- participant_scores %>%
  pivot_longer(cols = starts_with("Scoring"), names_to = "Scoring", values_to = "Score")

ggplot(long_data, aes(x = Score, y = ..density..)) +
  geom_histogram(binwidth = 10, alpha = 0.6) +
  geom_density(linewidth = 1) +
  facet_wrap(~ Scoring, scales = "free") +
  labs(title = "Histograms with Density Plots of Scoring Methods", x = "Score", y = "Density")

## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

make data separate

scoring_1_data <- long_format_data %>%
  filter(Scoring_Method == "Scoring_1") %>% 
  filter(Correct_Answer == "Yes")

scoring_2_data <- long_format_data %>%
  filter(Scoring_Method == "Scoring_2") %>%
  mutate(Score = as.ordered(Score))

scoring_3_data <- long_format_data %>%
  filter(Scoring_Method == "Scoring_3") %>%
  mutate(Score = as.ordered(Score))

scoring_4_data <- long_format_data %>%
  filter(Scoring_Method == "Scoring_4") %>%
  mutate(Score = as.ordered(Score))

scoring_5_data <- long_format_data %>%
  filter(Scoring_Method == "Scoring_5") %>%
  mutate(Score = as.ordered(Score))

datatable(
  scoring_1_data,
  filter='top', extensions = 'Scroller', class="compact",
  options = list(scrollY = 400, searching = T, paging = F,
                 columnDefs = list(list(className = 'dt-left', targets = "_all"))))

datatable(
  scoring_2_data,
  filter='top', extensions = 'Scroller', class="compact",
  options = list(scrollY = 400, searching = T, paging = F,
                 columnDefs = list(list(className = 'dt-left', targets = "_all"))))

## Warning in instance$preRenderHook(instance): It seems your data is too big for
## client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

datatable(
  scoring_3_data,
  filter='top', extensions = 'Scroller', class="compact",
  options = list(scrollY = 400, searching = T, paging = F,
                 columnDefs = list(list(className = 'dt-left', targets = "_all"))))

## Warning in instance$preRenderHook(instance): It seems your data is too big for
## client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

datatable(
  scoring_4_data,
  filter='top', extensions = 'Scroller', class="compact",
  options = list(scrollY = 400, searching = T, paging = F,
                 columnDefs = list(list(className = 'dt-left', targets = "_all"))))

## Warning in instance$preRenderHook(instance): It seems your data is too big for
## client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

datatable(
  scoring_5_data,
  filter='top', extensions = 'Scroller', class="compact",
  options = list(scrollY = 400, searching = T, paging = F,
                 columnDefs = list(list(className = 'dt-left', targets = "_all"))))

## Warning in instance$preRenderHook(instance): It seems your data is too big for
## client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

previous_scoring_data <- combined_data %>%
  mutate(
    Hit = ifelse(Response == "Yes" & Correct_Answer == "Yes", 1, 0),
    Miss = ifelse(Response == "No" & Correct_Answer == "Yes", 1, 0),
    False_Alarm = ifelse(Response == "Yes" & Correct_Answer == "No", 1, 0),
    Correct_Rejection = ifelse(Response == "No" & Correct_Answer == "No", 1, 0))

participant_rates <- previous_scoring_data %>%
  group_by(Participant) %>%
  summarise(
    h = sum(Hit) / sum(Correct_Answer == "Yes"),
    f = sum(False_Alarm) / sum(Correct_Answer == "No"))

participant_rates <- participant_rates %>%
  mutate(
    delta_m = (h - f) / (1 - f) - (f / h),
    I_SDT = 1 - (4 * h * (1 - f) - 2 * (h - f) * (1 + h - f)) / (4 * h * (1 - f) - (h - f) * (1 + h - f)))

describe(participant_rates) %>%
  kable(format = "html") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = F)

combined_data <- participant_rates %>%
  inner_join(participant_scores, by = "Participant")

cor_matrix <- combined_data %>%
  select(h, f, delta_m, I_SDT, Scoring_1, Scoring_2, Scoring_3, Scoring_4, Scoring_5) %>%
  cor()

sjPlot::tab_corr(
  cor_matrix, triangle = "lower", 
  title = "Correlation Matrix",
  var.labels = c("Hit Rate", "False Alarm Rate", "Delta_m", "I_SDT", "Scoring_1", "Scoring_2", "Scoring_3", "Scoring_4", "Scoring_5"))

ggpairs(
  combined_data, columns = 2:10, 
  title = "Correlation Matrix",
  upper = list(continuous = "cor"),
  lower = list(continuous = "smooth"))

Bayesian Estimation

prior <- set_prior("normal(0, 3)", class = "sd", group = "Participant") +
  set_prior("normal(0, 3)", class = "sd", group = "Item")

fit_scoring_1 <- brm(
  Score ~ 1 + (1 | Item) + (1 | Participant),
  data = scoring_1_data,
  family = bernoulli(link = "logit"),
  prior = prior,
  save_pars = save_pars(all = TRUE),
  seed = 1234)

## Compiling Stan program...

## Start sampling

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess

## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

fit_scoring_2 <- brm(
  Score ~ 1 + (1 | Item) + (1 | Participant),
  data = scoring_2_data,
  family = bernoulli(link = "logit"),
  prior = prior,
  save_pars = save_pars(all = TRUE),
  seed = 1234)

## Compiling Stan program...

## recompiling to avoid crashing R session

## Start sampling

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

Stan Code

make_stancode(fit_scoring_1)

## // generated with brms 2.21.0
## functions {
## }
## data {
##   int<lower=1> N;  // total number of observations
##   array[N] int Y;  // response variable
##   // data for group-level effects of ID 1
##   int<lower=1> N_1;  // number of grouping levels
##   int<lower=1> M_1;  // number of coefficients per level
##   array[N] int<lower=1> J_1;  // grouping indicator per observation
##   // group-level predictor values
##   vector[N] Z_1_1;
##   // data for group-level effects of ID 2
##   int<lower=1> N_2;  // number of grouping levels
##   int<lower=1> M_2;  // number of coefficients per level
##   array[N] int<lower=1> J_2;  // grouping indicator per observation
##   // group-level predictor values
##   vector[N] Z_2_1;
##   int prior_only;  // should the likelihood be ignored?
## }
## transformed data {
## }
## parameters {
##   real Intercept;  // temporary intercept for centered predictors
##   vector<lower=0>[M_1] sd_1;  // group-level standard deviations
##   array[M_1] vector[N_1] z_1;  // standardized group-level effects
##   vector<lower=0>[M_2] sd_2;  // group-level standard deviations
##   array[M_2] vector[N_2] z_2;  // standardized group-level effects
## }
## transformed parameters {
##   vector[N_1] r_1_1;  // actual group-level effects
##   vector[N_2] r_2_1;  // actual group-level effects
##   real lprior = 0;  // prior contributions to the log posterior
##   r_1_1 = (sd_1[1] * (z_1[1]));
##   r_2_1 = (sd_2[1] * (z_2[1]));
##   lprior += student_t_lpdf(Intercept | 3, 0, 2.5);
##   lprior += normal_lpdf(sd_1 | 0, 3)
##     - 1 * normal_lccdf(0 | 0, 3);
##   lprior += normal_lpdf(sd_2 | 0, 3)
##     - 1 * normal_lccdf(0 | 0, 3);
## }
## model {
##   // likelihood including constants
##   if (!prior_only) {
##     // initialize linear predictor term
##     vector[N] mu = rep_vector(0.0, N);
##     mu += Intercept;
##     for (n in 1:N) {
##       // add more terms to the linear predictor
##       mu[n] += r_1_1[J_1[n]] * Z_1_1[n] + r_2_1[J_2[n]] * Z_2_1[n];
##     }
##     target += bernoulli_logit_lpmf(Y | mu);
##   }
##   // priors including constants
##   target += lprior;
##   target += std_normal_lpdf(z_1[1]);
##   target += std_normal_lpdf(z_2[1]);
## }
## generated quantities {
##   // actual population-level intercept
##   real b_Intercept = Intercept;
## }

make_stancode(fit_scoring_2)

## // generated with brms 2.21.0
## functions {
## }
## data {
##   int<lower=1> N;  // total number of observations
##   array[N] int Y;  // response variable
##   // data for group-level effects of ID 1
##   int<lower=1> N_1;  // number of grouping levels
##   int<lower=1> M_1;  // number of coefficients per level
##   array[N] int<lower=1> J_1;  // grouping indicator per observation
##   // group-level predictor values
##   vector[N] Z_1_1;
##   // data for group-level effects of ID 2
##   int<lower=1> N_2;  // number of grouping levels
##   int<lower=1> M_2;  // number of coefficients per level
##   array[N] int<lower=1> J_2;  // grouping indicator per observation
##   // group-level predictor values
##   vector[N] Z_2_1;
##   int prior_only;  // should the likelihood be ignored?
## }
## transformed data {
## }
## parameters {
##   real Intercept;  // temporary intercept for centered predictors
##   vector<lower=0>[M_1] sd_1;  // group-level standard deviations
##   array[M_1] vector[N_1] z_1;  // standardized group-level effects
##   vector<lower=0>[M_2] sd_2;  // group-level standard deviations
##   array[M_2] vector[N_2] z_2;  // standardized group-level effects
## }
## transformed parameters {
##   vector[N_1] r_1_1;  // actual group-level effects
##   vector[N_2] r_2_1;  // actual group-level effects
##   real lprior = 0;  // prior contributions to the log posterior
##   r_1_1 = (sd_1[1] * (z_1[1]));
##   r_2_1 = (sd_2[1] * (z_2[1]));
##   lprior += student_t_lpdf(Intercept | 3, 0, 2.5);
##   lprior += normal_lpdf(sd_1 | 0, 3)
##     - 1 * normal_lccdf(0 | 0, 3);
##   lprior += normal_lpdf(sd_2 | 0, 3)
##     - 1 * normal_lccdf(0 | 0, 3);
## }
## model {
##   // likelihood including constants
##   if (!prior_only) {
##     // initialize linear predictor term
##     vector[N] mu = rep_vector(0.0, N);
##     mu += Intercept;
##     for (n in 1:N) {
##       // add more terms to the linear predictor
##       mu[n] += r_1_1[J_1[n]] * Z_1_1[n] + r_2_1[J_2[n]] * Z_2_1[n];
##     }
##     target += bernoulli_logit_lpmf(Y | mu);
##   }
##   // priors including constants
##   target += lprior;
##   target += std_normal_lpdf(z_1[1]);
##   target += std_normal_lpdf(z_2[1]);
## }
## generated quantities {
##   // actual population-level intercept
##   real b_Intercept = Intercept;
## }

Model Results

summary(fit_scoring_1)

##  Family: bernoulli 
##   Links: mu = logit 
## Formula: Score ~ 1 + (1 | Item) + (1 | Participant) 
##    Data: scoring_1_data (Number of observations: 10692) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Multilevel Hyperparameters:
## ~Item (Number of levels: 108) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.62      0.11     1.41     1.87 1.00      452     1013
## 
## ~Participant (Number of levels: 99) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.59      0.12     1.37     1.85 1.02      369      720
## 
## Regression Coefficients:
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.12      0.22    -0.58     0.31 1.02      157      217
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

summary(fit_scoring_2)

##  Family: bernoulli 
##   Links: mu = logit 
## Formula: Score ~ 1 + (1 | Item) + (1 | Participant) 
##    Data: scoring_2_data (Number of observations: 17820) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Multilevel Hyperparameters:
## ~Item (Number of levels: 180) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.94      0.11     1.75     2.17 1.02      286      696
## 
## ~Participant (Number of levels: 99) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.94      0.07     0.81     1.09 1.01      464     1011
## 
## Regression Coefficients:
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     1.15      0.17     0.83     1.49 1.02      117      395
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Fit Statistics

extract item and person parameters

ranef_scoring_1 <- ranef(fit_scoring_1)
ranef_scoring_2 <- ranef(fit_scoring_2)

item parameter

item_pars_scoring_1 <- ranef_scoring_1$Item
item_pars_scoring_2 <- ranef_scoring_2$Item

person parameter

person_pars_scoring_1 <- ranef_scoring_1$Participant
person_pars_scoring_2 <- ranef_scoring_2$Participant

coef_scoring_1 <- coef(fit_scoring_1)$Participant
coef_scoring_2 <- coef(fit_scoring_2)$Participant

# 項目の推定値を抽出
item_pars_scoring_1_estimate <- item_pars_scoring_1[, "Estimate", 1]

# 項目名を取得
item_names <- dimnames(item_pars_scoring_1)[[1]]

# データフレームに変換
item_data <- data.frame(Parameter = item_pars_scoring_1_estimate, Item = item_names, Type = "Item")

# 受験者の推定値を抽出
person_pars_scoring_1_estimate <- person_pars_scoring_1[, "Estimate", 1]

# 受験者名を取得
person_names <- dimnames(person_pars_scoring_1)[[1]]

# データフレームに変換
person_data <- data.frame(Parameter = person_pars_scoring_1_estimate, Participant = person_names, Type = "Person")

# 列名を統一
colnames(item_data)[colnames(item_data) == "Item"] <- "Name"
colnames(person_data)[colnames(person_data) == "Participant"] <- "Name"

# データを結合
all_data <- rbind(item_data, person_data)

# Extract the first dimension of item_pars_scoring_1, which corresponds to "Estimate"
item_data <- data.frame(Parameter = item_pars_scoring_1[, "Estimate", 1], 
                        Name = rownames(item_pars_scoring_1), 
                        Type = "Item")

# Extract the first dimension of person_pars_scoring_1, which corresponds to "Estimate"
person_data <- data.frame(Parameter = person_pars_scoring_1[, "Estimate", 1], 
                          Name = rownames(person_pars_scoring_1), 
                          Type = "Person")

# Combine the person and item data
combined_data <- rbind(item_data, person_data)

# ヒストグラムを左右に分けて描画
ggplot(combined_data, aes(x = Parameter, fill = Type)) +
  geom_histogram(data = subset(combined_data, Type == "Item"), aes(y = ..count..),
                 binwidth = 0.5, alpha = 0.7) +
  geom_histogram(data = subset(combined_data, Type == "Person"), aes(y = -..count..),
                 binwidth = 0.5, alpha = 0.7) +
  labs(x = "θ(Logit)", y = "Frequency") +
  coord_flip() +
  theme_bw()

## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

predict the probability of correct response

predicted_probs_1 <- predict(fit_scoring_1, type = "response")
predicted_probs_2 <- predict(fit_scoring_2, type = "response")

calculate residuals

numeric_score_1 <- as.numeric(as.character(fit_scoring_1$data$Score))
numeric_score_2 <- as.numeric(as.character(fit_scoring_2$data$Score))

residuals_1 <- numeric_score_1 - predicted_probs_1[, "Estimate"]
residuals_2 <- numeric_score_2 - predicted_probs_2[, "Estimate"]

item information function

info_function_1 <- predicted_probs_1[, "Estimate"] * (1 - predicted_probs_1[, "Estimate"])
info_function_2 <- predicted_probs_2[, "Estimate"] * (1 - predicted_probs_2[, "Estimate"])

outfit and infit statistics

\[ Outfit=\frac{1}{N}\sum_{i=1}^{N} \frac{(U_i - P_{ij})^2}{P_{ij}(1 - P_{ij})} \]

\[ Infit= \frac{\sum_{i=1}^{N} (U_{ij} - P_{ij})^2}{\sum_{i=1}^{N} P_{ij}(1 - P_{ij})} \]

calculate_fit_stats_1 <- function(residuals_1, info_function_1) {
  N <- length(residuals_1)
  
  outfit_ms <- mean(residuals_1^2 / info_function_1)
  outfit_sd <- sqrt((sum(1/info_function_1^2) / N^2) - 4/N)
  
  infit_ms <- sum(residuals_1^2) / sum(info_function_1)
  infit_sd <- sqrt((sum(info_function_1) - 4 * sum(info_function_1^2)) / sum(info_function_1)^2)
  
  outfit_t <- ((outfit_ms^(1/3) - 1) * (3/outfit_sd)) + (outfit_sd/3)
  infit_t <- ((infit_ms^(1/3) - 1) * (3/infit_sd)) + (infit_sd/3)
  
  return(c(outfit_ms = outfit_ms, outfit_t = outfit_t,
           infit_ms = infit_ms, infit_t = infit_t))
}

calculate_fit_stats_2 <- function(residuals_2, info_function_2) {
  N <- length(residuals_2)
  
  outfit_ms <- mean(residuals_2^2 / info_function_2)
  outfit_sd <- sqrt((sum(1/info_function_2^2) / N^2) - 4/N)
  
  infit_ms <- sum(residuals_2^2) / sum(info_function_2)
  infit_sd <- sqrt((sum(info_function_2) - 4 * sum(info_function_2^2)) / sum(info_function_2)^2)
  
  outfit_t <- ((outfit_ms^(1/3) - 1) * (3/outfit_sd)) + (outfit_sd/3)
  infit_t <- ((infit_ms^(1/3) - 1) * (3/infit_sd)) + (infit_sd/3)
  
  return(c(outfit_ms = outfit_ms, outfit_t = outfit_t,
           infit_ms = infit_ms, infit_t = infit_t))
}

item_fit_stats_1 <- fit_scoring_1$data %>%
  group_by(Item) %>%
  summarise(
    fit_stats_1 = list(calculate_fit_stats_1(residuals_1[cur_group_rows()], info_function_1[cur_group_rows()])),
    n = n()
  ) %>%
  tidyr::unnest_wider(fit_stats_1)

person_fit_stats_1 <- fit_scoring_1$data %>%
  group_by(Participant) %>%
  summarise(
    fit_stats_1 = list(calculate_fit_stats_1(residuals_1[cur_group_rows()], info_function_1[cur_group_rows()])),
    n = n()
  ) %>%
  tidyr::unnest_wider(fit_stats_1)

item_fit_stats_2 <- fit_scoring_2$data %>%
  group_by(Item) %>%
  summarise(
    fit_stats_2 = list(calculate_fit_stats_2(residuals_2[cur_group_rows()], info_function_2[cur_group_rows()])),
    n = n()
  ) %>%
  tidyr::unnest_wider(fit_stats_2)

person_fit_stats_2 <- fit_scoring_2$data %>%
  group_by(Participant) %>%
  summarise(
    fit_stats_2 = list(calculate_fit_stats_2(residuals_2[cur_group_rows()], info_function_2[cur_group_rows()])),
    n = n()
  ) %>%
  tidyr::unnest_wider(fit_stats_2)

ggplot(person_fit_stats_1, aes(x = infit_ms)) +
  geom_histogram(binwidth = 0.05, alpha = 0.7) +
  geom_vline(xintercept = c(0.5, 1.5), linetype = "dashed", color = "red") +
  labs(title = "Distribution of Infit Mean Square (Participants): Scoring_1",
       x = "Infit Mean Square", y = "Count")

ggplot(item_fit_stats_1, aes(x = infit_ms)) +
  geom_histogram(binwidth = 0.05, alpha = 0.7) +
  geom_vline(xintercept = c(0.5, 1.5), linetype = "dashed", color = "red") +
  labs(title = "Distribution of Infit Mean Square (Items): Scoring_1",
       x = "Infit Mean Square", y = "Count")

ggplot(item_fit_stats_1, aes(x = infit_ms, y = outfit_ms)) +
  geom_point(alpha = 0.7) +
  geom_text(aes(label = Item), vjust = -0.5, hjust = 0.5, size = 3) +
  geom_hline(yintercept = c(0.5, 1.5), linetype = "dashed", color = "red") +
  geom_vline(xintercept = c(0.5, 1.5), linetype = "dashed", color = "red") +
  labs(title = "Infit vs Outfit Mean Square (Items): Scoring_1",
       x = "Infit Mean Square", y = "Outfit Mean Square")

ggplot(person_fit_stats_2, aes(x = infit_ms)) +
  geom_histogram(binwidth = 0.05, alpha = 0.7) +
  geom_vline(xintercept = c(0.5, 1.5), linetype = "dashed", color = "red") +
  labs(title = "Distribution of Infit Mean Square (Participants): Scoring_2",
       x = "Infit Mean Square", y = "Count")

ggplot(item_fit_stats_2, aes(x = infit_ms)) +
  geom_histogram(binwidth = 0.05, alpha = 0.7) +
  geom_vline(xintercept = c(0.5, 1.5), linetype = "dashed", color = "red") +
  labs(title = "Distribution of Infit Mean Square (Items): Scoring_2",
       x = "Infit Mean Square", y = "Count")

ggplot(item_fit_stats_2, aes(x = infit_ms, y = outfit_ms)) +
  geom_point(alpha = 0.7) +
  geom_text(aes(label = Item), vjust = -0.5, hjust = 0.5, size = 3) +
  geom_hline(yintercept = c(0.5, 1.5), linetype = "dashed", color = "red") +
  geom_vline(xintercept = c(0.5, 1.5), linetype = "dashed", color = "red") +
  labs(title = "Infit vs Outfit Mean Square (Items): Scoring_2",
       x = "Infit Mean Square", y = "Outfit Mean Square")

Refit the model with removing outliers

list of items to remove

items_to_remove <- c("topic", "quotation", "invention", "discretion", "hignall", "descript")

New Scoring_1

scoring_1_data_filtered <- scoring_1_data %>%
  filter(!Item %in% items_to_remove)

(n_removed <- sum(scoring_1_data$Item %in% items_to_remove))

## [1] 297

(n_remaining <- nrow(scoring_1_data_filtered))

## [1] 10395

(n_unique_items_before <- length(unique(scoring_1_data$Item)))

## [1] 108

(n_unique_items_after <- length(unique(scoring_1_data_filtered$Item)))

## [1] 105

scoring_2_data_filtered <- scoring_2_data %>%
  filter(!Item %in% items_to_remove)

(n_removed <- sum(scoring_2_data$Item %in% items_to_remove))

## [1] 495

(n_remaining <- nrow(scoring_2_data_filtered))

## [1] 17325

(n_unique_items_before <- length(unique(scoring_2_data$Item)))

## [1] 180

(n_unique_items_after <- length(unique(scoring_2_data_filtered$Item)))

## [1] 175

fit_scoring_1 <- brm(
  Score ~ 1 + (1 | Item) + (1 | Participant),
  data = scoring_1_data_filtered,
  family = bernoulli(link = "logit"),
  prior = prior,
  save_pars = save_pars(all = TRUE),
  seed = 1234)

## Compiling Stan program...

## Trying to compile a simple C file

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
##          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1

## Start sampling

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000744 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 7.44 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 16.542 seconds (Warm-up)
## Chain 1:                14.488 seconds (Sampling)
## Chain 1:                31.03 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000447 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 4.47 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 14.184 seconds (Warm-up)
## Chain 2:                13.993 seconds (Sampling)
## Chain 2:                28.177 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000435 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 4.35 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 13.507 seconds (Warm-up)
## Chain 3:                13.859 seconds (Sampling)
## Chain 3:                27.366 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000443 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 4.43 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 14.203 seconds (Warm-up)
## Chain 4:                13.899 seconds (Sampling)
## Chain 4:                28.102 seconds (Total)
## Chain 4:

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess

## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

fit_scoring_2 <- brm(
  Score ~ 1 + (1 | Item) + (1 | Participant),
  data = scoring_2_data_filtered,
  family = bernoulli(link = "logit"),
  prior = prior,
  save_pars = save_pars(all = TRUE),
  seed = 1234)

## Compiling Stan program...
## Trying to compile a simple C file

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
##          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1

## Start sampling

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.001184 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 11.84 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 24.91 seconds (Warm-up)
## Chain 1:                22.579 seconds (Sampling)
## Chain 1:                47.489 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000703 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 7.03 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 24.22 seconds (Warm-up)
## Chain 2:                22.957 seconds (Sampling)
## Chain 2:                47.177 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000712 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 7.12 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 27.191 seconds (Warm-up)
## Chain 3:                23.319 seconds (Sampling)
## Chain 3:                50.51 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000716 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 7.16 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 24.279 seconds (Warm-up)
## Chain 4:                22.843 seconds (Sampling)
## Chain 4:                47.122 seconds (Total)
## Chain 4:

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

Stan Code

make_stancode(fit_scoring_1)

## // generated with brms 2.21.0
## functions {
## }
## data {
##   int<lower=1> N;  // total number of observations
##   array[N] int Y;  // response variable
##   // data for group-level effects of ID 1
##   int<lower=1> N_1;  // number of grouping levels
##   int<lower=1> M_1;  // number of coefficients per level
##   array[N] int<lower=1> J_1;  // grouping indicator per observation
##   // group-level predictor values
##   vector[N] Z_1_1;
##   // data for group-level effects of ID 2
##   int<lower=1> N_2;  // number of grouping levels
##   int<lower=1> M_2;  // number of coefficients per level
##   array[N] int<lower=1> J_2;  // grouping indicator per observation
##   // group-level predictor values
##   vector[N] Z_2_1;
##   int prior_only;  // should the likelihood be ignored?
## }
## transformed data {
## }
## parameters {
##   real Intercept;  // temporary intercept for centered predictors
##   vector<lower=0>[M_1] sd_1;  // group-level standard deviations
##   array[M_1] vector[N_1] z_1;  // standardized group-level effects
##   vector<lower=0>[M_2] sd_2;  // group-level standard deviations
##   array[M_2] vector[N_2] z_2;  // standardized group-level effects
## }
## transformed parameters {
##   vector[N_1] r_1_1;  // actual group-level effects
##   vector[N_2] r_2_1;  // actual group-level effects
##   real lprior = 0;  // prior contributions to the log posterior
##   r_1_1 = (sd_1[1] * (z_1[1]));
##   r_2_1 = (sd_2[1] * (z_2[1]));
##   lprior += student_t_lpdf(Intercept | 3, 0, 2.5);
##   lprior += normal_lpdf(sd_1 | 0, 3)
##     - 1 * normal_lccdf(0 | 0, 3);
##   lprior += normal_lpdf(sd_2 | 0, 3)
##     - 1 * normal_lccdf(0 | 0, 3);
## }
## model {
##   // likelihood including constants
##   if (!prior_only) {
##     // initialize linear predictor term
##     vector[N] mu = rep_vector(0.0, N);
##     mu += Intercept;
##     for (n in 1:N) {
##       // add more terms to the linear predictor
##       mu[n] += r_1_1[J_1[n]] * Z_1_1[n] + r_2_1[J_2[n]] * Z_2_1[n];
##     }
##     target += bernoulli_logit_lpmf(Y | mu);
##   }
##   // priors including constants
##   target += lprior;
##   target += std_normal_lpdf(z_1[1]);
##   target += std_normal_lpdf(z_2[1]);
## }
## generated quantities {
##   // actual population-level intercept
##   real b_Intercept = Intercept;
## }

make_stancode(fit_scoring_2)

## // generated with brms 2.21.0
## functions {
## }
## data {
##   int<lower=1> N;  // total number of observations
##   array[N] int Y;  // response variable
##   // data for group-level effects of ID 1
##   int<lower=1> N_1;  // number of grouping levels
##   int<lower=1> M_1;  // number of coefficients per level
##   array[N] int<lower=1> J_1;  // grouping indicator per observation
##   // group-level predictor values
##   vector[N] Z_1_1;
##   // data for group-level effects of ID 2
##   int<lower=1> N_2;  // number of grouping levels
##   int<lower=1> M_2;  // number of coefficients per level
##   array[N] int<lower=1> J_2;  // grouping indicator per observation
##   // group-level predictor values
##   vector[N] Z_2_1;
##   int prior_only;  // should the likelihood be ignored?
## }
## transformed data {
## }
## parameters {
##   real Intercept;  // temporary intercept for centered predictors
##   vector<lower=0>[M_1] sd_1;  // group-level standard deviations
##   array[M_1] vector[N_1] z_1;  // standardized group-level effects
##   vector<lower=0>[M_2] sd_2;  // group-level standard deviations
##   array[M_2] vector[N_2] z_2;  // standardized group-level effects
## }
## transformed parameters {
##   vector[N_1] r_1_1;  // actual group-level effects
##   vector[N_2] r_2_1;  // actual group-level effects
##   real lprior = 0;  // prior contributions to the log posterior
##   r_1_1 = (sd_1[1] * (z_1[1]));
##   r_2_1 = (sd_2[1] * (z_2[1]));
##   lprior += student_t_lpdf(Intercept | 3, 0, 2.5);
##   lprior += normal_lpdf(sd_1 | 0, 3)
##     - 1 * normal_lccdf(0 | 0, 3);
##   lprior += normal_lpdf(sd_2 | 0, 3)
##     - 1 * normal_lccdf(0 | 0, 3);
## }
## model {
##   // likelihood including constants
##   if (!prior_only) {
##     // initialize linear predictor term
##     vector[N] mu = rep_vector(0.0, N);
##     mu += Intercept;
##     for (n in 1:N) {
##       // add more terms to the linear predictor
##       mu[n] += r_1_1[J_1[n]] * Z_1_1[n] + r_2_1[J_2[n]] * Z_2_1[n];
##     }
##     target += bernoulli_logit_lpmf(Y | mu);
##   }
##   // priors including constants
##   target += lprior;
##   target += std_normal_lpdf(z_1[1]);
##   target += std_normal_lpdf(z_2[1]);
## }
## generated quantities {
##   // actual population-level intercept
##   real b_Intercept = Intercept;
## }

Model Results

summary(fit_scoring_1)

##  Family: bernoulli 
##   Links: mu = logit 
## Formula: Score ~ 1 + (1 | Item) + (1 | Participant) 
##    Data: scoring_1_data_filtered (Number of observations: 10395) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Multilevel Hyperparameters:
## ~Item (Number of levels: 105) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.56      0.12     1.35     1.80 1.02      312      892
## 
## ~Participant (Number of levels: 99) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.62      0.12     1.39     1.87 1.00      377      753
## 
## Regression Coefficients:
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.14      0.22    -0.59     0.27 1.02      155      350
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

summary(fit_scoring_2)

##  Family: bernoulli 
##   Links: mu = logit 
## Formula: Score ~ 1 + (1 | Item) + (1 | Participant) 
##    Data: scoring_2_data_filtered (Number of observations: 17325) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Multilevel Hyperparameters:
## ~Item (Number of levels: 175) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.92      0.10     1.72     2.13 1.01      399      800
## 
## ~Participant (Number of levels: 99) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.97      0.08     0.83     1.13 1.01      497     1071
## 
## Regression Coefficients:
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     1.13      0.18     0.78     1.46 1.01      128      272
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Model Comparison

The Leave-One-Out Cross-Validation (LOO-CV) method involves training a model on a dataset excluding one data point, then making a prediction on the excluded data point. This process is repeated for each data point in the dataset.

Consider a dataset \[ D = \{(x₁, y₁), \dots, (xₙ, yₙ)\} \]).

The LOO-CV estimate is expressed as follows:

\[ \text{LOO-CV} = -2 \sum_{i=1}^n \log p(y_i | x_i, D_{(-i)}) \]

loo_cv_scoring_1 <- brms::loo(fit_scoring_1)
loo_cv_scoring_2 <- brms::loo(fit_scoring_2)

print(loo_cv_scoring_1, simplify = FALSE)

## 
## Computed from 4000 by 10395 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -4675.5  55.4
## p_loo       192.9   2.9
## looic      9350.9 110.9
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [1.1, 2.9]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

print(loo_cv_scoring_2, simplify = FALSE)

## 
## Computed from 4000 by 17325 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -7154.5  79.8
## p_loo       256.8   4.4
## looic     14309.1 159.6
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [1.1, 2.9]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

The Widely Applicable Information Criterion (WAIC) is one of the information criteria used to evaluate the predictive performance of Bayesian models. It is asymptotically equivalent to Leave-One-Out Cross-Validation (LOO-CV) and is widely used due to its relatively easy computation.

WAIC is defined by the following formula:

\[ \text{WAIC} = -2(\text{lppd} - p_{\text{WAIC}}) \]

• lppd is the log pointwise predictive density.

• p_WAIC is the effective number of parameters (a penalty term for model complexity).

waic_1 <- brms::waic(fit_scoring_1)

## Warning: 
## 1 (0.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.

waic_2 <- brms::waic(fit_scoring_2)

## Warning: 
## 14 (0.1%) p_waic estimates greater than 0.4. We recommend trying loo instead.

print(waic_1$estimates)

##             Estimate         SE
## elpd_waic -4675.0887  55.427852
## p_waic      192.5204   2.882035
## waic       9350.1773 110.855704

print(waic_2$estimates)

##             Estimate        SE
## elpd_waic -7153.8132  79.78581
## p_waic      256.0621   4.35038
## waic      14307.6264 159.57162

Visualization

plot(fit_scoring_1, combo = c("trace", "dens_overlay"), ask = F)

plot(fit_scoring_2, combo = c("trace", "dens_overlay"), ask = F)

pp_check(
  object = fit_scoring_1,
  ndraws = 1000,
  type = "stat",
  stat = "mean") +
  theme_test() +
  ylab("Frequency")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(
  object = fit_scoring_1,
  ndraws = 100) +
  theme_test() +
  ylab("Density")

pp_check(
  object = fit_scoring_2,
  ndraws = 1000,
  type = "stat",
  stat = "mean") +
  theme_test() +
  ylab("Frequency")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(
  object = fit_scoring_2,
  ndraws = 100) +
  theme_test() +
  ylab("Density")

ranef_scoring_1$Participant[, , "Intercept"] %>%
  as_tibble() %>%
  rownames_to_column() %>%
  arrange(Estimate) %>%
  mutate(rowname = seq_len(n())) %>%
  ggplot(aes(rowname, Estimate, ymin = Q2.5, ymax = Q97.5)) +
  geom_pointrange(alpha = 0.7) +
  coord_flip() +
  labs(x = "Person Number (Sorted)")

ranef_scoring_1$Item[, , "Intercept"] %>%
  as_tibble() %>%
  rownames_to_column() %>%
  arrange(Estimate) %>%
  mutate(rowname = seq_len(n())) %>%
  ggplot(aes(rowname, Estimate, ymin = Q2.5, ymax = Q97.5)) +
  geom_pointrange(alpha = 0.7) +
  coord_flip() +
  labs(x = "Item Number (Sorted)")

ranef_scoring_2$Participant[, , "Intercept"] %>%
  as_tibble() %>%
  rownames_to_column() %>%
  arrange(Estimate) %>%
  mutate(rowname = seq_len(n())) %>%
  ggplot(aes(rowname, Estimate, ymin = Q2.5, ymax = Q97.5)) +
  geom_pointrange(alpha = 0.7) +
  coord_flip() +
  labs(x = "Person Number (Sorted)")

ranef_scoring_2$Item[, , "Intercept"] %>%
  as_tibble() %>%
  rownames_to_column() %>%
  arrange(Estimate) %>%
  mutate(rowname = seq_len(n())) %>%
  ggplot(aes(rowname, Estimate, ymin = Q2.5, ymax = Q97.5)) +
  geom_pointrange(alpha = 0.7) +
  coord_flip() +
  labs(x = "Item Number (Sorted)")

calculate_p_correct <- function(theta, b, a = 1) {
  1 / (1 + exp(-a * (theta + b)))
}

calculate_item_info <- function(theta, b, a = 1) {
  p <- calculate_p_correct(theta, b, a)
  a^2 * p * (1 + p)
}

calculate_test_info <- function(theta, item_params) {
  sapply(theta, function(t) sum(mapply(calculate_item_info, t, item_params)))
}

create_plot_data <- function(theta, selected_items, ranef_data_1, ranef_data_2) {
  item_diff_1 <- ranef_data_1$Item[selected_items, , 1, drop = TRUE]
  item_diff_2 <- ranef_data_2$Item[selected_items, , 1, drop = TRUE]
  
  plot_data <- expand.grid(theta = theta, item = selected_items, scoring = c("Scoring_1", "Scoring_2"))
  
  plot_data <- plot_data %>%
    mutate(
      difficulty = case_when(
        scoring == "Scoring_1" ~ item_diff_1[match(as.character(item), selected_items), "Estimate"],
        scoring == "Scoring_2" ~ item_diff_2[match(as.character(item), selected_items), "Estimate"]),
      lower_ci = case_when(
        scoring == "Scoring_1" ~ item_diff_1[match(as.character(item), selected_items), "Q2.5"],
        scoring == "Scoring_2" ~ item_diff_2[match(as.character(item), selected_items), "Q2.5"]),
      upper_ci = case_when(
        scoring == "Scoring_1" ~ item_diff_1[match(as.character(item), selected_items), "Q97.5"],
        scoring == "Scoring_2" ~ item_diff_2[match(as.character(item), selected_items), "Q97.5"]),
      p = mapply(calculate_p_correct, theta, difficulty),
      p_lower = mapply(calculate_p_correct, theta, lower_ci),
      p_upper = mapply(calculate_p_correct, theta, upper_ci),
      info = mapply(calculate_item_info, theta, difficulty),
      info_lower = mapply(calculate_item_info, theta, lower_ci),
      info_upper = mapply(calculate_item_info, theta, upper_ci))
  
  return(plot_data)
}

create_tif_data <- function(theta, ranef_data_1, ranef_data_2) {
  calculate_tif_with_ci <- function(theta, ranef_data) {
    item_params <- ranef_data$Item[, c("Estimate", "Q2.5", "Q97.5"), 1]
    test_info <- apply(item_params, 2, function(params) calculate_test_info(theta, params))
    data.frame(theta = theta, 
               test_info = test_info[, "Estimate"],
               lower_ci = test_info[, "Q2.5"],
               upper_ci = test_info[, "Q97.5"])
  }
  
  bind_rows(
    calculate_tif_with_ci(theta, ranef_data_1) %>% mutate(scoring = "Scoring_1"),
    calculate_tif_with_ci(theta, ranef_data_2) %>% mutate(scoring = "Scoring_2"))
}

create_diff_data <- function(ranef_data_1, ranef_data_2) {
  bind_rows(
    data.frame(ranef_data_1$Item[, , 1], item = rownames(ranef_data_1$Item), scoring = "Scoring_1"),
    data.frame(ranef_data_2$Item[, , 1], item = rownames(ranef_data_2$Item), scoring = "Scoring_2")) %>%
    select(item, Estimate, Q2.5, Q97.5, scoring)
}

create_comparison_plots <- function(ranef_scoring_1, ranef_scoring_2, num_items = 5) {
  common_items <- intersect(rownames(ranef_scoring_1$Item), rownames(ranef_scoring_2$Item))
  selected_items <- sample(common_items, min(num_items, length(common_items)))
  
  theta <- seq(-6, 6, length.out = 100)
  
  plot_data <- create_plot_data(theta, selected_items, ranef_scoring_1, ranef_scoring_2)
  tif_data <- create_tif_data(theta, ranef_scoring_1, ranef_scoring_2)
  diff_data <- create_diff_data(ranef_scoring_1, ranef_scoring_2)
  
  icc_plot <- create_icc_plot(plot_data)
  iic_plot <- create_iic_plot(plot_data)
  tif_plot <- create_tif_plot(tif_data)
  diff_plot <- create_diff_plot(diff_data)
  
  list(icc = icc_plot, iic = iic_plot, tif = tif_plot, difficulty = diff_plot)
}

create_icc_plot <- function(plot_data) {
  ggplot(plot_data, aes(x = theta, y = p, color = item, linetype = scoring)) +
    geom_line() +
    geom_ribbon(aes(ymin = p_lower, ymax = p_upper, fill = item), alpha = 0.1, color = NA) +
    labs(title = "Item Characteristic Curves",
         x = "Latent Trait (θ)", y = "Probability of Correct Response") 
}

create_iic_plot <- function(plot_data) {
  ggplot(plot_data, aes(x = theta, y = info, color = item, linetype = scoring)) +
    geom_line() +
    geom_ribbon(aes(ymin = info_lower, ymax = info_upper, fill = item), alpha = 0.1, color = NA) +
    labs(title = "Item Information Curves",
         x = "Latent Trait (θ)", y = "Item Information")
}

create_tif_data <- function(theta, ranef_data_1, ranef_data_2) {
  calculate_tif_with_ci <- function(theta, ranef_data) {
    item_params <- ranef_data$Item[, c("Estimate", "Q2.5", "Q97.5"), 1, drop = FALSE]
    test_info <- apply(item_params, 2, function(params) calculate_test_info(theta, params))
    data.frame(theta = theta, 
               test_info = test_info[, "Estimate"],
               lower_ci = test_info[, "Q2.5"],
               upper_ci = test_info[, "Q97.5"])
  }
  
  bind_rows(
    calculate_tif_with_ci(theta, ranef_data_1) %>% mutate(scoring = "Scoring_1"),
    calculate_tif_with_ci(theta, ranef_data_2) %>% mutate(scoring = "Scoring_2"))
}

create_tif_plot <- function(tif_data) {
  ggplot(tif_data, aes(x = theta, y = test_info, color = scoring)) +
    geom_line() +
    geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci, fill = scoring), alpha = 0.1, color = NA) +
    labs(title = "Test Information Function",
         x = "Latent Trait (θ)", y = "Test Information")
}

create_diff_data <- function(ranef_data_1, ranef_data_2) {
  bind_rows(
    as.data.frame(ranef_data_1$Item[, , 1]) %>% 
      mutate(item = rownames(ranef_data_1$Item), scoring = "Scoring_1"),
    as.data.frame(ranef_data_2$Item[, , 1]) %>% 
      mutate(item = rownames(ranef_data_2$Item), scoring = "Scoring_2")) %>%
    select(item, Estimate, Q2.5, Q97.5, scoring)
}

create_diff_plot <- function(diff_data) {
  ggplot(diff_data, aes(x = Estimate, y = reorder(item, Estimate), color = scoring)) +
    geom_point(position = position_dodge(width = 0.5)) +
    geom_errorbarh(aes(xmin = Q2.5, xmax = Q97.5), height = 0.2, position = position_dodge(width = 0.5)) +
    labs(title = "Item Difficulties Comparison",
         x = "Difficulty (b-parameter)", y = "Item")
}

create_comparison_plots <- function(ranef_scoring_1, ranef_scoring_2, num_items = 5) {
  common_items <- intersect(rownames(ranef_scoring_1$Item), rownames(ranef_scoring_2$Item))
  selected_items <- sample(common_items, min(num_items, length(common_items)))
  
  theta <- seq(-6, 6, length.out = 100)
  
  plot_data <- create_plot_data(theta, selected_items, ranef_scoring_1, ranef_scoring_2)
  tif_data <- create_tif_data(theta, ranef_scoring_1, ranef_scoring_2)
  diff_data <- create_diff_data(ranef_scoring_1, ranef_scoring_2)
  
  icc_plot <- create_icc_plot(plot_data)
  iic_plot <- create_iic_plot(plot_data)
  tif_plot <- create_tif_plot(tif_data)
  diff_plot <- create_diff_plot(diff_data)
  
  list(icc = icc_plot, iic = iic_plot, tif = tif_plot, difficulty = diff_plot)
}

create_comparison_plots(ranef_scoring_1, ranef_scoring_2)

## $icc

## 
## $iic

## 
## $tif

## 
## $difficulty

theta_range <- seq(-6, 6, length.out = 100)

item_info <- function(b, theta) {
  p_theta <- 1 / (1 + exp(-(theta + b)))
  return(p_theta * (1 - p_theta))
}

# Calculate item information for each item
calc_item_info <- function(b_participant, theta_range) {
  info_list <- lapply(seq_along(theta_range), function(t) {
    theta = theta_range[t]
    mean_info = sum(sapply(b_participant[, "Estimate"], item_info, theta = theta))
    lower_info = sum(sapply(b_participant[, "Q2.5"], item_info, theta = theta))
    upper_info = sum(sapply(b_participant[, "Q97.5"], item_info, theta = theta))
    return(c(mean = mean_info, lower = lower_info, upper = upper_info))
  })
  return(do.call(rbind, info_list))
}

# Extract participant ability parameters for each scoring method
b_participant_scoring_1 <- coef(fit_scoring_1)$Participant[, ,"Intercept"]
b_participant_scoring_2 <- coef(fit_scoring_2)$Participant[, ,"Intercept"]

# Calculate item information for each scoring method
info_results_scoring_1 <- calc_item_info(b_participant_scoring_1, theta_range)
info_results_scoring_2 <- calc_item_info(b_participant_scoring_2, theta_range)

# Convert the results to data frames for plotting
info_df_scoring_1 <- as.data.frame(info_results_scoring_1)
info_df_scoring_1$theta <- theta_range

info_df_scoring_2 <- as.data.frame(info_results_scoring_2)
info_df_scoring_2$theta <- theta_range

info_df_scoring_1 <- info_df_scoring_1 %>%
  mutate(lower = pmin(lower, mean),
         upper = pmax(upper, mean))

info_df_scoring_2 <- info_df_scoring_2 %>%
  mutate(lower = pmin(lower, mean),
         upper = pmax(upper, mean))

# maximum test information and theta
max_scoring_1 <- max(info_df_scoring_1$upper)
max_scoring_2 <- max(info_df_scoring_2$upper)

max_scoring_1_row <- info_df_scoring_1[info_df_scoring_1$upper == max_scoring_1, ]
max_scoring_2_row <- info_df_scoring_2[info_df_scoring_2$upper == max_scoring_2, ]

max_theta_scoring_1 <- max_scoring_1_row$theta
max_theta_scoring_2 <- max_scoring_2_row$theta

max_mean_scoring_1 <- max_scoring_1_row$mean
max_mean_scoring_2 <- max_scoring_2_row$mean

cat("Theta for maximum test information in Scoring 1:", max_theta_scoring_1, "\nMean for maximum test information in Scoring 1:", max_mean_scoring_1)

## Theta for maximum test information in Scoring 1: -0.7878788 
## Mean for maximum test information in Scoring 1: 16.63914

cat("\nTheta for maximum test information in Scoring 2:", max_theta_scoring_2, "\nMean for maximum test information in Scoring 2:", max_mean_scoring_2)

## 
## Theta for maximum test information in Scoring 2: -1.151515 
## Mean for maximum test information in Scoring 2: 20.59561

# Add a scoring method identifier to each data frame
info_df_scoring_1 <- info_df_scoring_1 %>% mutate(Scoring_Method = "Scoring_1")
info_df_scoring_2 <- info_df_scoring_2 %>% mutate(Scoring_Method = "Scoring_2")

combined_info_df <- bind_rows(info_df_scoring_1, info_df_scoring_2)

ggplot(combined_info_df, aes(x = theta, y = mean, color = Scoring_Method, fill = Scoring_Method)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2) +
  labs(title = "Item Information Function by Scoring Method", x = "Theta", y = "Information")

Sensitivity Analysis

prior distributions’ SDs changes from 1 to 5

prior_settings <- list(
  prior_normal_1 = set_prior("normal(0, 1)", class = "sd", group = "Participant") +
    set_prior("normal(0, 1)", class = "sd", group = "Item"),
  prior_normal_2 = set_prior("normal(0, 2)", class = "sd", group = "Participant") +
    set_prior("normal(0, 2)", class = "sd", group = "Item"),
  prior_normal_3 = set_prior("normal(0, 3)", class = "sd", group = "Participant") +
    set_prior("normal(0, 3)", class = "sd", group = "Item"),
  prior_normal_4 = set_prior("normal(0, 4)", class = "sd", group = "Participant") +
    set_prior("normal(0, 4)", class = "sd", group = "Item"),
  prior_normal_5 = set_prior("normal(0, 5)", class = "sd", group = "Participant") +
    set_prior("normal(0, 5)", class = "sd", group = "Item"))

scoring data list

scoring_data_list <- list(
  scoring_1 = scoring_1_data,
  scoring_2 = scoring_2_data)

model_results <- list()
ranef_results <- list()

for (scoring_method in names(scoring_data_list)) {
  for (prior_setting in names(prior_settings)) {
    if (scoring_method %in% c("scoring_1", "scoring_2")) {
      family <- bernoulli(link = "logit")
    } else {
      family <- cumulative(link = "logit")
    }
    
    model_1pl <- brm(
      Score ~ 1 + (1 | Item) + (1 | Participant),
      data = scoring_data_list[[scoring_method]],
      family = family,
      prior = prior_settings[[prior_setting]],
      seed = 1234,
      control = list(adapt_delta = 0.95)
    )

    model_results[[scoring_method]][[prior_setting]] <- model_1pl
    ranef_results[[scoring_method]][[prior_setting]] <- ranef(model_1pl)
  }
}

## Compiling Stan program...

## Trying to compile a simple C file

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
##          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1

## Start sampling

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000822 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 8.22 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 20.153 seconds (Warm-up)
## Chain 1:                14.831 seconds (Sampling)
## Chain 1:                34.984 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000439 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 4.39 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 19.74 seconds (Warm-up)
## Chain 2:                14.257 seconds (Sampling)
## Chain 2:                33.997 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000501 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 5.01 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 23.492 seconds (Warm-up)
## Chain 3:                14.287 seconds (Sampling)
## Chain 3:                37.779 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000439 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 4.39 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 18.861 seconds (Warm-up)
## Chain 4:                14.356 seconds (Sampling)
## Chain 4:                33.217 seconds (Total)
## Chain 4:

## Warning: The largest R-hat is 1.06, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess

## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

## Compiling Stan program...
## Trying to compile a simple C file

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
##          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1

## Start sampling

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000683 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 6.83 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 23.132 seconds (Warm-up)
## Chain 1:                14.537 seconds (Sampling)
## Chain 1:                37.669 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000446 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 4.46 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 22.082 seconds (Warm-up)
## Chain 2:                16.44 seconds (Sampling)
## Chain 2:                38.522 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000448 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 4.48 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 25.114 seconds (Warm-up)
## Chain 3:                14.52 seconds (Sampling)
## Chain 3:                39.634 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000445 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 4.45 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 20.593 seconds (Warm-up)
## Chain 4:                15.332 seconds (Sampling)
## Chain 4:                35.925 seconds (Total)
## Chain 4:

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

## Compiling Stan program...
## Trying to compile a simple C file

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
##          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1

## Start sampling

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000668 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 6.68 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 21.691 seconds (Warm-up)
## Chain 1:                15.182 seconds (Sampling)
## Chain 1:                36.873 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000457 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 4.57 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 20.28 seconds (Warm-up)
## Chain 2:                29.518 seconds (Sampling)
## Chain 2:                49.798 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000448 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 4.48 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 21.079 seconds (Warm-up)
## Chain 3:                16.807 seconds (Sampling)
## Chain 3:                37.886 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000467 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 4.67 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 21.695 seconds (Warm-up)
## Chain 4:                14.73 seconds (Sampling)
## Chain 4:                36.425 seconds (Total)
## Chain 4:

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

## Compiling Stan program...
## Trying to compile a simple C file

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
##          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1

## Start sampling

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.00087 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 8.7 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 21.046 seconds (Warm-up)
## Chain 1:                14.215 seconds (Sampling)
## Chain 1:                35.261 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000451 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 4.51 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 19.122 seconds (Warm-up)
## Chain 2:                19.521 seconds (Sampling)
## Chain 2:                38.643 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.00045 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 4.5 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 20.459 seconds (Warm-up)
## Chain 3:                14.852 seconds (Sampling)
## Chain 3:                35.311 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000443 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 4.43 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 19.795 seconds (Warm-up)
## Chain 4:                17.106 seconds (Sampling)
## Chain 4:                36.901 seconds (Total)
## Chain 4:

## Warning: The largest R-hat is 1.05, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess

## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

## Compiling Stan program...
## Trying to compile a simple C file

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
##          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1

## Start sampling

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000694 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 6.94 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 22.139 seconds (Warm-up)
## Chain 1:                20.705 seconds (Sampling)
## Chain 1:                42.844 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000442 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 4.42 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 20.279 seconds (Warm-up)
## Chain 2:                17.982 seconds (Sampling)
## Chain 2:                38.261 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.00046 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 4.6 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 24.231 seconds (Warm-up)
## Chain 3:                14.12 seconds (Sampling)
## Chain 3:                38.351 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000438 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 4.38 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 18.399 seconds (Warm-up)
## Chain 4:                22.054 seconds (Sampling)
## Chain 4:                40.453 seconds (Total)
## Chain 4:

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

## Compiling Stan program...
## Trying to compile a simple C file

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
##          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1

## Start sampling

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.001318 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 13.18 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 41.988 seconds (Warm-up)
## Chain 1:                29.4 seconds (Sampling)
## Chain 1:                71.388 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000744 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 7.44 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 30.884 seconds (Warm-up)
## Chain 2:                23.876 seconds (Sampling)
## Chain 2:                54.76 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000711 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 7.11 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 48.712 seconds (Warm-up)
## Chain 3:                24.082 seconds (Sampling)
## Chain 3:                72.794 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000816 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 8.16 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 33.668 seconds (Warm-up)
## Chain 4:                24.684 seconds (Sampling)
## Chain 4:                58.352 seconds (Total)
## Chain 4:

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

## Compiling Stan program...
## Trying to compile a simple C file

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
##          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1

## Start sampling

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.001294 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 12.94 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 38.341 seconds (Warm-up)
## Chain 1:                23.639 seconds (Sampling)
## Chain 1:                61.98 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000719 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 7.19 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 40.139 seconds (Warm-up)
## Chain 2:                23.856 seconds (Sampling)
## Chain 2:                63.995 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000743 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 7.43 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 37.568 seconds (Warm-up)
## Chain 3:                23.827 seconds (Sampling)
## Chain 3:                61.395 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000989 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 9.89 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 36.009 seconds (Warm-up)
## Chain 4:                48.098 seconds (Sampling)
## Chain 4:                84.107 seconds (Total)
## Chain 4:

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

## Compiling Stan program...
## Trying to compile a simple C file

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
##          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1

## Start sampling

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.001215 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 12.15 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 37.601 seconds (Warm-up)
## Chain 1:                46.219 seconds (Sampling)
## Chain 1:                83.82 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000715 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 7.15 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 34.111 seconds (Warm-up)
## Chain 2:                40.922 seconds (Sampling)
## Chain 2:                75.033 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000743 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 7.43 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 42.314 seconds (Warm-up)
## Chain 3:                28.261 seconds (Sampling)
## Chain 3:                70.575 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000713 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 7.13 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 40.093 seconds (Warm-up)
## Chain 4:                23.136 seconds (Sampling)
## Chain 4:                63.229 seconds (Total)
## Chain 4:

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

## Compiling Stan program...
## Trying to compile a simple C file

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
##          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1

## Start sampling

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.00118 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 11.8 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 35.888 seconds (Warm-up)
## Chain 1:                23.12 seconds (Sampling)
## Chain 1:                59.008 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000715 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 7.15 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 33.025 seconds (Warm-up)
## Chain 2:                46.117 seconds (Sampling)
## Chain 2:                79.142 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.00074 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 7.4 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 35.563 seconds (Warm-up)
## Chain 3:                23.113 seconds (Sampling)
## Chain 3:                58.676 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000726 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 7.26 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 1074.23 seconds (Warm-up)
## Chain 4:                25.552 seconds (Sampling)
## Chain 4:                1099.79 seconds (Total)
## Chain 4:

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

## Compiling Stan program...
## Trying to compile a simple C file

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
##          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1

## Start sampling

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.001335 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 13.35 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 35.652 seconds (Warm-up)
## Chain 1:                22.983 seconds (Sampling)
## Chain 1:                58.635 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000721 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 7.21 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 394.248 seconds (Warm-up)
## Chain 2:                23.104 seconds (Sampling)
## Chain 2:                417.352 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000706 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 7.06 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 43.946 seconds (Warm-up)
## Chain 3:                45.759 seconds (Sampling)
## Chain 3:                89.705 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000715 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 7.15 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 36.885 seconds (Warm-up)
## Chain 4:                165.951 seconds (Sampling)
## Chain 4:                202.836 seconds (Total)
## Chain 4:

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

extract_and_combine <- function(ranef_list) {
  do.call(rbind, lapply(names(ranef_list), function(prior_setting) {
    do.call(rbind, lapply(names(ranef_list[[prior_setting]]), function(group) {
      df <- as.data.frame(ranef_list[[prior_setting]][[group]])
      df$prior_setting <- prior_setting
      df$group <- group
      df$item <- rownames(df)
      df
    }))
  }))
}

combined_results <- lapply(ranef_results, extract_and_combine)

combined_data <- do.call(rbind, lapply(names(combined_results), function(scoring_method) {
  df <- combined_results[[scoring_method]]
  df$scoring_method <- scoring_method
  df
}))

ggplot(combined_data, aes(x = prior_setting, y = Estimate.Intercept, group = interaction(item, scoring_method), color = scoring_method)) +
  geom_line(alpha = 0.5) +
  geom_point(size = 1) +
  labs(title = "Item Difficulties under Different Priors and Scoring Methods",
       x = "Prior Setting",
       y = "Estimate (Intercept)") +
  theme(legend.title = element_blank(),
        legend.position = "bottom")

ggplot(combined_data, aes(x = prior_setting, y = Estimate.Intercept, group = item, color = item)) +
  geom_line(alpha = 0.7) +
  geom_point(size = 1) +
  labs(title = "Item Difficulties under Different Priors and Scoring Methods",
       x = "Prior Setting",
       y = "Estimate (Intercept)") +
  theme(legend.position = "none") +
  facet_wrap(~ scoring_method, ncol = 1)

ggplot(combined_data, aes(x = item, y = Estimate.Intercept, color = scoring_method)) +
  geom_point(size = 2) +
  labs(title = "Item Difficulties under Different Priors and Scoring Methods",
       x = "Item",
       y = "Estimate (Intercept)") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1),
        legend.title = element_blank(),
        legend.position = "bottom") +
  facet_wrap(~ prior_setting, ncol = 2)

ggplot(combined_data, aes(x = item, y = scoring_method, fill = Estimate.Intercept)) +
  geom_tile() +
  labs(title = "Item Difficulties under Different Priors and Scoring Methods",
       x = "Item",
       y = "Scoring Method",
       fill = "Estimate (Intercept)") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1),
        legend.position = "bottom") +
  facet_wrap(~ prior_setting, ncol = 2)

ggplot(combined_data, aes(x = prior_setting, y = Estimate.Intercept, fill = prior_setting)) +
  geom_violin(alpha = 0.7) +
  geom_boxplot(width = 0.1, fill = "white", outlier.shape = NA) +
  labs(title = "Comparison of Item Difficulty Estimates under Different Priors",
       x = "Prior Setting",
       y = "Estimate (Intercept)") +
  theme(legend.position = "none") +
  facet_wrap(~ scoring_method, ncol = 1)

ggplot(combined_data, aes(x = Estimate.Intercept, y = prior_setting, fill = prior_setting)) +
  geom_density_ridges(alpha = 0.7) +
  labs(title = "Comparison of Item Difficulty Estimates under Different Priors",
       x = "Estimate (Intercept)",
       y = "Prior Setting") +
  theme(legend.position = "none") +
  facet_wrap(~ scoring_method, ncol = 1)

## Picking joint bandwidth of 0.483

## Picking joint bandwidth of 0.471

Multidimensional IRT

assume that hit and correct rejection are completely different

fit_scoring_2_mirt <- brm(
    Score ~ 1 + (1 | Item) + (0 + Correct_Rejection | Participant),
    data = scoring_2_data,
    family = bernoulli(link = "logit"),
    prior = prior,
    seed = 1234)

## Compiling Stan program...

## Trying to compile a simple C file

## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## #include <cmath>
##          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1

## Start sampling

## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.001218 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 12.18 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 67.538 seconds (Warm-up)
## Chain 1:                83.222 seconds (Sampling)
## Chain 1:                150.76 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.000779 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 7.79 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 65.899 seconds (Warm-up)
## Chain 2:                182.846 seconds (Sampling)
## Chain 2:                248.745 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0.000775 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 7.75 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 68.403 seconds (Warm-up)
## Chain 3:                49.551 seconds (Sampling)
## Chain 3:                117.954 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0.000785 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 7.85 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 84.493 seconds (Warm-up)
## Chain 4:                49.637 seconds (Sampling)
## Chain 4:                134.13 seconds (Total)
## Chain 4:

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess

## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

summary(fit_scoring_2_mirt)

##  Family: bernoulli 
##   Links: mu = logit 
## Formula: Score ~ 1 + (1 | Item) + (0 + Correct_Rejection | Participant) 
##    Data: scoring_2_data (Number of observations: 17820) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Multilevel Hyperparameters:
## ~Item (Number of levels: 180) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.82      0.14     1.57     2.13 1.01      420      922
## 
## ~Participant (Number of levels: 99) 
##                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Correct_Rejection)    19.99      1.67    16.94    23.38 1.00     3841
##                       Tail_ESS
## sd(Correct_Rejection)     3138
## 
## Regression Coefficients:
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -1.14      0.17    -1.48    -0.83 1.02      140      277
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

loo_cv_scoring_2_mirt <- brms::loo(fit_scoring_2_mirt)

## Warning: Found 2451 observations with a pareto_k > 0.7 in model
## 'fit_scoring_2_mirt'. We recommend to set 'moment_match = TRUE' in order to
## perform moment matching for problematic observations.

print(loo_cv_scoring_2_mirt, simplify = FALSE)

## 
## Computed from 4000 by 17820 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -6358.9  57.9
## p_loo       112.4   1.6
## looic     12717.7 115.8
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.4]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     15369 86.2%   1646    
##    (0.7, 1]   (bad)       2417 13.6%   <NA>    
##    (1, Inf)   (very bad)    34  0.2%   <NA>    
## See help('pareto-k-diagnostic') for details.

plot(fit_scoring_2_mirt, combo = c("trace", "dens_overlay"), ask = F)

pp_check(
  object = fit_scoring_2_mirt,
  ndraws = 1000,
  type = "stat",
  stat = "mean") +
  theme_test() +
  ylab("Frequency")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(
  object = fit_scoring_2_mirt,
  ndraws = 100) +
  theme_test() +
  ylab("Density")

ranef_scoring_2_mirt <- ranef(fit_scoring_2_mirt)

item_pars_scoring_2_mirt <- ranef_scoring_2_mirt$Item

person_pars_scoring_2_mirt <- ranef_scoring_2_mirt$Participant

coef_scoring_2_mirt <- coef(fit_scoring_2_mirt)$Participant

ranef_scoring_2_mirt$Participant[, , "Correct_Rejection"] %>%
  as_tibble() %>%
  rownames_to_column() %>%
  arrange(Estimate) %>%
  mutate(rowname = seq_len(n())) %>%
  ggplot(aes(rowname, Estimate, ymin = Q2.5, ymax = Q97.5)) +
  geom_pointrange(alpha = 0.7) +
  coord_flip() +
  labs(x = "Person Number (Sorted)")

ranef_scoring_2_mirt$Item[, , "Intercept"] %>%
  as_tibble() %>%
  rownames_to_column() %>%
  arrange(Estimate) %>%
  mutate(rowname = seq_len(n())) %>%
  ggplot(aes(rowname, Estimate, ymin = Q2.5, ymax = Q97.5)) +
  geom_pointrange(alpha = 0.7) +
  coord_flip() +
  labs(x = "Item Number (Sorted)")

theta_range <- seq(-6, 6, length.out = 100)

item_info <- function(b, theta) {
  p_theta <- 1 / (1 + exp(-(theta - b)))
  return(p_theta * (1 - p_theta))
}

calc_item_info <- function(b_participant, theta_range) {
  info_list <- lapply(seq_along(theta_range), function(t) {
    theta = theta_range[t]
    mean_info = sum(sapply(b_participant[, "Estimate"], item_info, theta = theta))
    lower_info = sum(sapply(b_participant[, "Q2.5"], item_info, theta = theta))
    upper_info = sum(sapply(b_participant[, "Q97.5"], item_info, theta = theta))
    return(c(mean = mean_info, lower = lower_info, upper = upper_info))
  })
  return(do.call(rbind, info_list))
}

# Extract participant ability parameters for each scoring method
b_participant_scoring_2_mirt <- coef(fit_scoring_2_mirt)$Participant[, ,"Correct_Rejection"]

# Calculate item information for each scoring method
info_results_scoring_2_mirt <- calc_item_info(b_participant_scoring_2_mirt, theta_range)

# Convert the results to data frames for plotting
info_df_scoring_2_mirt <- as.data.frame(info_results_scoring_2_mirt)
info_df_scoring_2_mirt$theta <- theta_range

info_df_scoring_2_mirt <- info_df_scoring_2_mirt %>%
  mutate(lower = pmin(lower, mean),
         upper = pmax(upper, mean))

max_scoring_2_mirt <- max(info_df_scoring_2_mirt$upper)
max_scoring_2_mirt_row <- info_df_scoring_2_mirt[info_df_scoring_2_mirt$upper == max_scoring_2_mirt, ]
max_theta_scoring_2_mirt <- max_scoring_2_mirt_row$theta
max_mean_scoring_2_mirt <- max_scoring_2_mirt_row$mean

cat("Theta for maximum test information in Scoring 2 with MIRT:", max_theta_scoring_2_mirt, "\nMean for maximum test information in Scoring 2 with MIRT:", max_mean_scoring_2_mirt)

## Theta for maximum test information in Scoring 2 with MIRT: 6 
## Mean for maximum test information in Scoring 2 with MIRT: 2.021271e-05

ggplot(info_df_scoring_2_mirt, aes(x = theta, y = mean)) +
  geom_line(color = "green") +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, fill = "green") +
  labs(title = "Item Information Function (Scoring 2 with MIRT)", x = "Theta", y = "Information")

# Add a scoring method identifier to each data frame
info_df_scoring_2_mirt <- info_df_scoring_2_mirt %>% mutate(Scoring_Method = "Scoring_2_mirt")

# Combine all data frames into one
combined_info_df_mirt <- bind_rows(info_df_scoring_2, info_df_scoring_2_mirt)

ggplot(combined_info_df_mirt, aes(x = theta, y = mean, color = Scoring_Method, fill = Scoring_Method)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2) +
  labs(title = "Item Information Function by Scoring Method", x = "Theta", y = "Information") +
  scale_color_manual(values = c("Scoring_2" = "red", "Scoring_2_mirt" = "green")) +
  scale_fill_manual(values = c("Scoring_2" = "red", "Scoring_2_mirt" = "green"))